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Abstract 

Through molecular dynamics, we study the d — 2,3 classical model of N cou- 
pled rotators (inertial XY model) assuming a coupling constant which decays 
with distance as r^" (a > 0). The total energy < H > is asymptotically oc NN 
with N = [N-^~°'^'^ — {a/d)]/[l — a/d], hence the model is thermodynamically 
extensive if a/d > 1 and nonextensive otherwise. We numerically show that, for 
energies above some threshold, the (appropriately scaled) maximum Lyapunov 
exponent is oc N^"^ where k is an universal (one and the same for c? = 1,2 
and 3, and all energies) function of a/d, which monotonically decreases from 
1/3 to zero when a/d increases from to 1, and identically vanishes above 1. 
These features are consistent with the nonextensive statistical mechanics sce- 
nario, where thermodynamic extensivity is associated with exponential mixing 
in phase space, whereas weaker (possibly power-law in the present case) mixing 
emerges at the N ~^ oo limit whenever nonextcnsivity is observed. 
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In the last few years there has been a noticeable interest in the study of 
the thermodynamics and statistical mechanics of anomalous systems exhibiting 
nonextensivity [1, 2]. A system is extensive if its energy and entropy, as func- 
tions of intensive internal parameters (e.g. temperature), grow linearly with the 
size of the system (i.e. with N, the number of its microscopic components). 
Nonextensivity can be brought into scene by long-range interactions. Prom 
the static point of view, nonextensivity means that thermodynamic quantities 
like the internal or the free energies per particle are not constant with N, but 
rather diverge in the thermodynamic limit N —>■ oo. Prom the dynamical point 
of view, one might think that long-range interactions induce persistent spatial 
and temporal correlations leading to the breaking of standard mixing and er- 
godicity properties, hence to a possible violation of the usual Boltzmann-Gibbs 
(BG) statistics. Classical Hamiltonian systems with many degrees of freedom 
are paradigmatic in the discussion of the foundations of equilibrium statistical 
mechanics. The maximal Lyapunov exponent (MLE) is a well established indi- 
cator of chaos. If it is positive, the system will gencrically be strongly chaotic 
and will satisfy the standard ergodic hypothesis attained through exponentially 
quick mixing. If it is instead zero non ergodicity could emerge, giving origin to 
less-than-exponcntial (typically power-law) mixing, hence to anomalous thcrmo- 
statistical behavior. It is then interesting to study maximal Lyapunov exponents 
in systems with long-range interactions, where nonextensive thermodynamical 
and dynamical anomalies arc expected. In a previous work [3] the maximal 
Lyapunov exponent of a one-dimensional system of N planar rotators coupled 
with interactions decaying as the inverse power a of their distances (see model 
(1) below) has been studied as a function of A^. For total energy above some 
threshold, the maximal Lyapunov exponent is proportional to where K{a) 

is a function which goes from k(0) > to zero while a increases from zero to 1 
(long range forces) and remains zero for a > 1 (short range forces) . In this letter 
we want to investigate if the exponent k, describing the weakening of the mixing 
properties of the dynamics in the long-range-interacting regime, is universal. In 
particular, we check the hypothesis that K{a,d) is a universal function of a/d, 
where d is the dimensionality of space. 

Let us consider a simple model of N planar rotators (XY spins) placed at 
the sites i = 1 ... TV of a square and a cubic lattice (i.e., dimension d = 2 and 3). 
The rotators have unit moment of inertia, angular momentum Lj and position 
specified by the angle 6i e [0, 27r]. Rotations occur inside an arbitrary reference 
plane. The model is described by the following classical hamiltonian in the 
conjugate canonical pairs {9i,Li}: 



where r^j is the distance between lattice site i and j with periodic boundary 
conditions in the form of nearest image convention; a > is a variable parameter 
specifying the range of the interaction, which is defined long-range if a < d and 
short-range if a > c?. 

This system is nonextensive in the case of long-range interactions [4]. To 
verify that, let us consider as an example a quantity related to the internal 
energy per rotator at T = (i. e., when K = and all rotators are equally 
oriented), namely the integral drr'^~^r~°' , which converges for a > d and 
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diverges for < a < d. In the latter case, the thermodynamics of the model 
is not well defined in the usual sense. However, it is possible to define energy 
quantities that are bounded in the thermodynamic limit if they are properly 
scaled (see [2, 5, 6]). It is interesting to remark that, in system (1) for < 
a < d, not only equilibrium (static) quantities, but also a quantity related to 
the local instability of the dynamics such as the maximal Lyapunov exponent 
might diverge with N [7] , and needs to be properly rcscalcd by a function of N 
if one wants to evaluate the degree of chaoticity as a function of the size of the 
system. 

An alternative procedure is to force the system to be extensive via a largely 
used mean- field-like rescaling [8] of the potential energy in (1) which then be- 
comes: 



The rescaling factor N in front of the potential energy V of hamiltonian (2), 
is a function of a, d, N and of the geometry of the lattice and has to be chosen 
to guarantee the existence of a bounded energy density U =< H'/N > in the 
thermodynamic limit. 

The model (1) has been introduced in [3] on a one-dimensional lattice and 
its thermodynamics has been studied by means of simulations in [9]. When 
a = and N — N the model is also called the Hamiltonian Mean-Field XY 
Model (HMF) [10], whose dynamics and thermodynamics have been extensively 
studied [11, 12, 13]. In this case the underlying lattice has no meaning and 
any dependence from its dimension d disappears. It is found that the HMF 
model has a phase transition from a ferromagnetic state to a paramagnetic one 
at the critical temperature = 0.5 corresponding to a critical energy density 
Uc = 0.75 (fcs = 1). In this letter we use [14] 



whose properties are shown in figure 1. (Unity has been added to the integral in 

the definition of A'^ in order to have the convenient limit A'^ ^ 1 when a/d oo.) 
Using a similar rescaling factor it was shown in [15] that model (2) and HMF 
model have the same thermodynamics, and in particular the same critical energy 
density. 

Chaoticity properties and maximal Lyapunov exponent (MLE) of HMF 
model have been studied by Latora, Rapisarda and Ruffo [13, 16, 17, 18]. They 
calculated the curve MLE versus U numerically for various A^ and analytically 
for A^ — *• oo; MLE goes to both when U tends to and when U tends to oo. 
In these two limits the system becomes integrable, being represented by a set 
of harmonic oscillators around the ferromagnetic fundamental state of aligned 
rotators in the former case and by free rotators in the latter. Finite-size ef- 
fects are present for U > Uc where MLE vanishes in the thermodynamic limit 
as ^ N^'^/'^. For fixed U < Uc the curve is smooth and presents a positive 
maximum at the critical energy Uc, where therefore in the thermodynamic limit 
exists a finite discontinuity. [19] 

The behavior of MLE in the one-dimensional model (2) has been studied, as 
said above, in [3], as a function of A^ and a at a fixed high energy (above Uc). It 
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has been shown that it goes to in the thermodynamic hmit when a < d, like 
in the HMF model, but it remains constant when a > d. The way the exponent 
vanishes is still of the form ~ N~'^ where now k is a function of a. The curve 
K{a) starts from the value 1/3 at a = and goes to when a 1~ . 

In the present work we extend this study to cases d > 1. We still expect 
that when a < d and U > Uc MLE vanishes in the thermodynamic limit with 
a power law of the type N~'^ with k ~ K{a, d) satisfying K{a = 0, d) = 1/3 [19] 
and vanishing when a — > cJ~. Two questions seem quite relevant. Is this curve 
K{a,d) universal in the sense that it only depends ona/d, or it also depends on 
the dimensionality d of the system ? Moreover, is this curve independent from 
the energy density U at which it is computed? In what follows we will positively 
answer to both questions. 

We have simulated the constant energy dynamics of the nonextensive hamil- 
tonian (1) and computed the maximal Lyapunov exponents. 

Hamiltonians H (eq. (1)) and H' (eq. (2)) are strictly connected. If we 
multiply H' by N, the potential energy V becomes V. Then, if we rescale time 
as t' = t^^N, and using the fact that L,'s involve first order time derivatives, we 
obtain that also the kinetic energy K' becomes K. Finally we get H = NH' . 

Integrating the equations of motion of H is identical to integrating those 
of NH' but for times longer by a factor \/iV- AH thermodynamic quantities 
that can be obtained as temporal averages of time dependent observables are 
not influenced by the rescaling of time. In particular, once one has T, V, U from 
H simulations, one can get the corresponding values for H' simply dividing by 
N . The evaluation of maximal Lyapunov exponents \'max ^f hamiltonian H' 
can also be obtained through rescaling of the one measured for H, noted Xmax 
hereafter. 

In order to calculate maximal Lyapunov exponents one must consider the 
limit, [20] 

Amax = lim 7 In = hm \{t) (4) 



with d{t) = \JY^i{50iY + {5LiY being the metric distance calculated from in- 
finitesimal displacements at time t obtained in turn through the double integra- 
tion of both the normal equations of motion 

Oi = Li (5) 



and the linearized ones 



SOi = 5L, (7) 
SL, = (8) 

(we report the ones derived from H, the ones for H' are exactly the same with 
the substitution rg- I'fjN). Because of formula (4) which involves 1/t, time 

rescaUng yields A^„j, = Xmax/V^. 

We computed \„iax for typical values of ~ 30 to 4000 and a fixed energy 
density U = 5.0 above the critical one Uc = 0.75. We made simulations on a 
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square and simple cubic lattices (i.e., d = 2,3) with unitary lattice step, for 
various values of the a/d ratio. 

Simulations started from an high temperature initial state: Oi's are randomly 
extracted from a imiform distribution in [0, 27r], momenta i^'s are extracted 
from a uniform distribution in [—0.5, 0.5] and then translated and rescaled to 
have zero total momentum and the desired total energy. We used the velocity- 
Verlet algorithm [21] with a time-step chosen to have a relative energy con- 
servation of 10"'* or better. Length of simulations were chosen looking at the 
asymptotic behavior of quantity X{t), eq. (4), where d{0) is randomly chosen 
(see [20]). We report in figure 2, the curves X'maxi^) various a in dimensions 
d = 2 and d = 3. 

These curves extend those of fig. 3 in [3]. In both cases we observe that, for 
growing N, if a > d then X'^^ax positive and constant, whereas if a < the 
maximal Lyapunov exponent tends to zero. 

We fitted the data with the following functional form: AJ^^^. oc N~'^; k being 
the slope in the log- log plots of fig. 2. 

In figure 3 we collect the slopes til cLS 3j function of the ratio a/d, for d = 1 
(data from [3]), d = 2 and d = 3. Remarkably, through the simple scaling a/d 
we show that the TV dependence of the rescaled maximal Lyapunov exponent 
of model (1) is universal. All curves start at a = from a value close to 
1/3, analytically and numerically found for the HMF model [19, 16]; when a/d 
increases from zero to unity all our data collapse into a single curve, and then 
remain zero for a/d greater than unity. All of our data are roughly fitted by 
the heuristic expression k= [1 — {a/d)^]/[3 + {a/d)'^/2]. 

Extrapolating for N going to infinity the results in fig. 2, it is interesting to 
note that the constant k behavior for a > d is consistent with the results found 
for the a = 00 (only first-neighboring coupling) model [22] . 

The universal function k shown in fig. 3 does not appear to depend on the 
energy density, provided this is greater than the critical one. For the case a = 
this has been proved analytically in [19] and numerically in [16]. We show in 
fig. 4 that this is also true for a. = 0.8 and d = 2, the asymptotic slope of A^^^, 
being sensibly the same for the two energy densities U = 5, U = 8. 

Summarizing our results, we can say that the classical hamiltonian H is non 
extensive for < a < d and extensive for a/d > 1. Through a N rescaling 
of its potential energy it can artificially be made extensive, thus merging with 
a long standing tradition in mean-field-like approaches. Above some critical 
value for the (rescaled) total energy per particle, the weakening of the mixing 
properties is universal; more precisely, the maximal Lyapunov exponents scale 
as N~'^, K being a function only of the ratio a/d (and not d nor U). In fact, 
it seems reasonable to expect that this K,{a/d) function is model- independent. 
More precisely, for total energies high enough, it is expected to have one and the 
same functional form K{a/d) for all classical models whose potential is nowhere 
singular and whose attractive tail decays as 1/r" . Consistently with this specu- 
lation, some specific gravitational and fluid models [23] do exhibit «;(0) = 1/3, in 
agreement with the present results for the XY models under analysis. An uni- 
versal K{a/d) fits in fact very well the generic expectation for such models within 
nonextensive statistical mechanics. Indeed, the available studies of hamiltonian 
systems including long range interactions provide a variety of indications that a 
crossover occurs, as a function of time, in the thermostatistical properties of the 
system. For times below a crossover time t{N), several properties are anoma- 
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lous (diffusion, velocity distribution and others), whereas they become normal 
for times much larger than r. The fact that t seems to diverge with N makes 
the anomalous regime to become the only observable one in the thermodynamic 
limit. It is precisely this regime which is focused on within that generalized 
formalism. For the specific case of the mixing properties, the present results 
suggest that a less-than-exponential mixing takes place until a crossover time 
after which the mixing becomes exponential, i.e., the usual one. This crossover 
time might well scale like l/X'max (or a positive power of it). Therefore, since 
^'max vanishes in the limit N ^ oo, we expect the crossover time to diverge, 
thus emerging the anomalous regime we were just discussing. 
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Figure 1: Scaling factor N(N, a/d): (a) versus a/d, for typical values of A*' (dotted 
line: the A/" — > oo asymptote); (b) versus A'', for typical values of a/d (dotted line: the 
a/d = 2 asymptote). Note the divergence of N in the thermodynamic limit {N — > oo), 
in the long-range case. 
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Figure 2: Curves X'max vs A'' (log-log scale) for d = 2 (upper part) and d = 3 
(lower part), U = 5. Solid lines are fits with the functional form (a — /{Ny. 
Consequently, X'^ax AT"'^'™''*' where K{a,d) = (1 — a/d)c{a,d) for < a < d and 
K(a, d) — for a > d; for a = d, X'^^^ is oxpoctod to vanish as a power of 1/ In A'". The 
ordinate indistinctively represents X'^^x corresponding to Hamiltonian (2) with U = 

H'/N = 5, or Xmax/V^ corresponding to Hamiltonian (1) with U = H/{NN) = 5. 
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Figure 3: The mixing weakening exponent n ys a/d for d= 1,2,3 {d — 1: from [3]; 
d = 2,3: present work). It describes the asymptotic A'^ behaviour of the maximal 
Lyapunov exponent X'max sit fixed energy above the critical one, i.e., X'max 1/N'^. 
The solid line is a guide to the eye consistent with universality. For a = we have 
[19] k(0) = 1/3 (Vd). 
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Figure 4: X',nax{N) for fixed a — 0.8, d = 2 and two different energies U. The 
asymptotic N behavior, for aU values of U above the critical one Uc{Q:,d) and < 
a/d < 1, appears to be X'max c< i/jv"'™''''-', where the proportionality coefficient 
decreases from a finite (a, d)-dependent positive value to zero when U increases from 
Uc{a, d) to infinity. 
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